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Abstract 


A nonoverlapping domain decomposition technique applied to a finite 
difference method is presented for the numerical solution of the forward- 
backward heat equation in the case of one-dimension. While the previous at- 
tempts in dealing with this problem have been based on an iterative domain 
decomposition scheme, the current work avoids iterations. Also a physical 
matching condition is suggested to avoid difficulties caused by the interface 
boundary nodes. Furthermore, we obtain a square system of equations. In 
addition, the convergence and stability of the proposed method are investi- 
gated. Some numerical experiments are given to show the effectiveness of the 
proposed method. 
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1 Introduction 


The analysis and numerical solution of the forward-backward heat equa- 
tion (FBHE) has been under consideration for more than three decades; see 
[1, 8, 12, 20]. This problem appears in many applications such as bound- 
ary layer problems in fluid dynamics and steady state computation, plasma 
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physics, and studies of the propagation of an electron beam through the 
solar corona (see [7] and the references therein). It also arises in computa- 
tional fluid dynamics and randomly accelerated particle problem, and some 
more applications are mentioned in [1]. The forward-backward problem has 
been solved by various methods such as finite difference method (FDM) 
[7, 6, 13, 23], transformation to a system of first order differential equa- 
tions [17], least square approach [3], Galerkin finite element [2, 10, 9, 16], 
and radial basis function meshless method [4]. 


The application of the domain decomposition method (DDM) together 
with FDM for this problem in the case of one-dimension is widely seen in 
the literature (see, for example, [7, 22]). This method was also applied to 
the two-dimensional FBHE, and convergence of the iterative method was 
considered in [23]. 


In addition, the theoretical aspects of the underlying problem have been 
studied by some authors [5, 14, 19]. Daoud [7] analyzed the convergence of 
the overlapping Schwarz waveform relaxation method for solving the FBHE. 


Kuznetsov [15] proved the uniqueness of the entropy solution for a nonlin- 
ear forward-backward parabolic equation with Dirichlet boundary condition 
and showed that initial and final conditions must be formulated in the form 
of inequalities. Paronetto [20] considered the existence and uniqueness of the 
solution for elliptic problems with a small parameter and proved that the so- 
lution converges to the solution of a mixed type equation, elliptic-parabolic, 
parabolic both forward and backward. In the above-mentioned methods, 
the domain is usually divided into two subdomains, each of which is asso- 
ciated with either a standard forward or backward problem. A numerical 
method is then applied to each subproblem. A primary approximate solution 
is assumed on the interior boundary followed by solving each subproblem 
and updating the interface boundary solution iteratively. While the previous 
DDM-FDM has been applied based on an iterative scheme, the current work 
is based on applying a nonoverlapping DDM and employing an FDM without 
using iterations. Some advantages of the proposed method read as follows. 
Firstly, an initial approximate solution on the interface boundary, which may 
affect the convergence, is not required. Secondly, the iteration, which may 
considerably increase the computational costs, is avoided. In the suggested 
method, the domain is divided into two nonoverlapping subdomains, and a 
forward or a backward finite difference formulation is applied to each sub- 
problem according to the sort of the initial conditions. The local algebraic 
equations produced by the subproblems are then assembled in a global system 
via employing a physical matching condition on the interface and obtaining 
a square system of equations by removing the virtual boundary nodes. This 
paper is organized as follows: The forward-backward equation is introduced 
in Section 2. The DDM will be applied to the underlying problem in Section 
3. The convergence and stability of the proposed method will be discussed 
in Section 4. Finally some numerical results will be presented in Section 5. 
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2 The forward-backward equation 


A boundary value problem of a forward-backward parabolic equation in one- 
dimension is introduced as 


LUt — Ure = f(x,t), (x,t) € Q = (1,1) x (0,1), 

), for all ¢ € [0,1], 

), for all ¢ € [0,1], 

u(a,0) = uo(a), for all x € [0,1], 

u(a,1) =ur(x), for all # € [-1,0], (1) 


where uo(a), ui(x), g_i(t), and gi(t) are known functions with uo(1) = gi (0) 
and wu (—1) = g_i(1). 


The existence, uniqueness, and stability of the above forward-backward 
problem have been considered by some authors [1, 2, 3, 18]. Authors of [1] 
used the energy method to prove the existence and uniqueness of a weak so- 
lution to FBHE on some special spaces equipped with some suitable norms. 
Lu and Wen [18] studied the existence and uniqueness of a weak solution 
to (1) on a certain Hilbert space. Aziz and Liu [3] reduced (1) to a first- 
order symmetric-positive system of differential equations and proved exis- 
tence and uniqueness to the solution of the equivalent problem using the 
theory of symmetric-positive systems previously presented by Friedrichs [11] 
for a strong solution, under certain smoothness assumptions. 


In order to apply the FDM, we produce a set of grids as follows: Let 
1 1 
h=—, let x; = ih fori=0,+1,..., (M —1), let 7 = 5, and let t; = jr 


M 
for 7 =0,...,N. 


3 Domain decomposition method 


In this section, the DDM is applied to the forward backward problem in two 
steps. First, the domain is divided into two subdomains and the algebraic 
equations are formed for each subregion. Then, the local formulations are 
assembled to obtain the global system of equations. 
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Figure 1: The domain partitioning and the pattern of the finite difference schemes. 


3.1 Domain partitioning 


Let the domain 2 be divided into two nonoverlapping subregions, 0; = 
(—1,0) x (0,1) and Q2 = (0,1) x (0,1) with the real boundaries T, and T, 
and an interface boundary I’; in between (see Figure 1). In order to deal 
with this problem, we apply a forward and a backward difference scheme 
to the domains Q; and Qe, respectively. Employing the first order forward 
formulation for time derivative and the second order central difference scheme 
for the spatial derivative at the mesh points in 01, we obtain 


jp ies 7 Ui — Wit 7 Qui j + Ui-1,7 


rT _ he ; + fig, 
0<j<N-1, -M4+1<i<-1, 
U-m,j =9-107), OSG<N. (2) 


where w;,; denotes the approximation of u(ih,j7) and fi,; = f(ih,jr). 

Using the first order backward formulation and the second order central 
difference scheme for the mesh points in Q2, the following equations are 
achieved: 


Fy ec — Uig — UWi+1,7+1 — Qui j+1 + Uj-1,j41 


. h2 + Jigits 


O<j<N-1, 1<i<M-1, 
umj=MT), OSISN, 
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The pattern of the finite difference schemes used in (2) and (3) is shown in 
Figure 1. 

Note that the above equations produce two separate underdetermined 
systems of equations, as the boundary conditions on the interface [7 are not 
prescribed. In the next subsection, we propose a method in which the local 
equations are assembled in a global system with a square coefficient matrix. 


3.2 Global formulation 


We now globalize the problem by assembling the local equations constructed 
by the subproblems. As mentioned above, the systems obtained by equations 
(2) and (3) are underdetermined as the number of unknowns is equal to 
2|N(M — 2) +N — 1], whereas there are only 2N(M — 2) equations, taking 
into account the boundary conditions imposed on the interface boundary 
have different values for two subproblems. 

Let uf®) , k = 1,2, be the solutions of the subproblems associated with OQ, 
at (0,77), 7 =1,...,.N. In order to achieve a system with a square coefficient 
matrix, we use a sort of physical matching condition on the interface, that 
is, Uo,j = uf!) = Thee which guarantees a continuous solution on the virtual 
boundary. We now make a square system by eliminating the nodal values 
located on the interface. To do so, we utilize the central difference scheme 
for spatial derivative at the mesh points on I’; in order to make interaction 


between two subproblems, that is, 


u1,j; — 2U0,j + U-1,5 . 
(Se EN) = fy, 1S iSN-1, (4) 


Note that the first term of the governing equation in (1) vanishes since on 
[';, « = 0 holds. By combining this equation with (2) and (3) in the case of 
7 = 1 andi = -—1, we can eliminate the nodal values on the line x = 0 from 
the eventual system of equations. 

To this end, we consider the first equation of (3) in the case of ¢ = 1, that 
is, 


(2) 
Ujt1 — Uj — U2ajt1 — 2U1j41 + Up j41 
h— = oS he + fijg, OS FS N-1, (5) 


where the superscript indicates the subdomain considered. 
Replacing j by j +1 in (4) and solving it for uo j41 = TOL we obtain 


h? 1 1 : 
Uy = Glos t suas t stig, OS7SN—1. (6) 


Substituting (6) into (5), we find that 
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3 1 
jo = “13 = bie _ aa il os { a . (7) 
T h2 Pg ts 
$0 esd om hy, 
Repeating a similar process for the line i = —1 in the first equation of (2), 


we achieve a similar equation without presence of the term uo,;4+1 as follows: 


3 
U—-2,j — FU-1,j + 5U1,3 
U-1,j41 — U-1,3 a a eee (8) 
h s = i = he + 5 fo,35 QS gs ag Ne 


Assembling equations (2) and (3) with the revised equations (7) and (8), 
a linear system of size 2N(M — 2) is obtained as PU = F, where P = 


Pry Pir Ur Fr 
Us , f= ‘ 


A, Ar 
B B 
Ba ; Pag |" , 
Ay '. AR 
By, Ay Br Ap 
Cr, 0 Cr 0 
Pur = . : Pat = oe ie 
L 0 Cr 0 
0 0 
2r —(-M + 1)h —r 
—r 2r —(-M 4+ 2)h —-r 
AL = , 
—r 2r —(—2)h —r 
—r 3r — (—h) 
2r—(-—M+1)h —r 
—r 2r —(-M + 2)h —-r 
Ap = % ? 
—r 2r —(—h) 
0 (—M +1)h 
os By = ; 
—£ —h 
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3r +h -r 
—r wr+2h—r 
Ar = ; é ‘ ’ 
—r 2r+(M—2)h —r 
—r 2r+(M—1)h 
2r+h -—-r 
—r 2r+2h—r 
Apr al é F A ’ 
—r 2r+(M—2)h —r 
—r 2r+(M-—1)h 
-£ _h 
Cr = “iess ; Br = nai ) 
0 (-M + 1)h, 


z 
where r = on 


The indices ZL and R denote the portions of the matrices associated with 
the left and right subregions, respectively. For instance, Pr represents the 
portions associated with coupling the nodal values of subproblem 1 (left side) 
to the nodal values of subproblem 2 (right side). 

Having solved the above linear system and obtained the interior solution 
of the subproblems, the interface boundary solution can be found by using 


(6). 


4 Convergence 


In order to show the convergence of the proposed method, we use Lax’s equiv- 
alence theorem [21]. To do so, we need to seek the stability condition and 
consistency of this method. We first check the stability. 


Theorem 1. The finite difference equations (2)—(4) are unconditionally sta- 
ble. 


Proof. Substituting up, = e°?"E% into the difference equation (2), we find 
that 


ph(é — 1) = r€(2.cos(Bh) — 2) = r€(—4r sin?()) 


Since ph < 1, we can write 


e(ar sin?) +1)>1. 
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Using the stability condition, |€| < 1, we have 


1 
—_, Bh... <|é| <1: 
Ar sin oe. +1 


Clearly 0 < € < 1 for all r > O and all 6. Therefore, the equations are 
unconditionally stable. For equation (3), the same result can be obtained. 


The consistency of the method is considered in the following theorem. 


Theorem 2. The finite difference equations (2)—(4) are consistent and they 
have the local truncation error of O(hk) + O(h?). 


Proof. Let Fi;(u) = 0 represent the difference equation approximating the 
partial differential equation (1) at the (i, 7)th mesh point, with exact solution 
u, and let U be the exact solution of equation (1). Then we have the local 
truncation error T;; at the mesh point (th, j7) as follows: 


Ti; = Fi;(U) 
oe Og g34 — Ug Opi ga — 20 ja + Gea ght = 
= jp J ch = thst fli,7 +1). 
By Taylor’s expansion, we have 
_, OU 0?U _ 
Ti = th( aL eget (5 5 Negeri — 70,9 +4) 


{hk GU kU h? AU WOU ’ 
‘o ee °3 OB 12 Ox! 20 dx : 


Hence 
Tij = O(hk) + O(h*). 


4.1 Higher order finite difference schemes 


The accuracy of the solution can be improved by employing higher order 
finite difference approximations. To do so, we first apply the previous finite 
difference formulation to the nodes located on the lines x = th for i= M—-1 
and i = —M +1 as follows: 


jplist Tg _ tg T 2uij + Wi-1,5 


T h2 + fas: 


Now there is sufficient number of nodal values to use the following scheme 
for the other nodes: 
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Uijt2 — 4Uij41 + 3ui,j 
27 
_ SUi42,5 1 16u;—-1,; — 80u;,; + 16u;;-1 — Ui j—2 
a 12h2 


ih 


+ fixe (9) 


Applying this formula improves the accuracy of the solution, taking into 
account that the order of the scheme (9) is O(hk?) + O(h*). As will be 
observed, the numerical errors are remarkably declined in the case of using 
this formula. 


5 Numerical results 


We now present some numerical results in order to show the effectiveness of 
the proposed method. 

Consider equation (1) with uo(x) = 0, ui(a%) = 0, g_i(t) = 0, gi(t) = 0, 
and 


2x (x? — 1)t[(t — 1)? — 4a? + t(t — 1)] — 2t?[(t — 1)? — 24x? + 4], 
_ x>0, te [0,1], 
PO) = 9 on(a? — 1)(¢ — 1)(22 — t — 4a?) — 2(t — 1)2(2 — 242? +4), 
x<0, t€ [0,1], 


The exact solution of the above problem is given by 


(x? — 1)t?[(t — 1)? — 427], x >0,t€ [0,1], 
u(x,t) = 


(x? — 1)(t? -t —4a7)(E-1)?,. x <0,t€ [0,1]. 


To measure the accuracy of the results, the maximum error (Max error) and 
the root-mean-square error (RMSE) are employed, that is, 


N 
So (tj a uj)?, 
j=l 


Z| 


Na 
Maz error = max| tj —u;|, RMSE= 
I= 


where N is the number of nodes and u; and t; denote the exact and approxi- 
mate solutions at the jth node. The numerical solution of the above equation 
was obtained in two cases: Scheme (I), the FDM proposed in (2) and (3), 
and scheme (II), the FDM suggested in (9). The numerical errors are listed 
in Tables 1 and 2, respectively, for the cases (I) and (II). The agreement of 
the approximate and the exact solutions together with the configuration of 
the error functions is shown in Figures 2 and 3, respectively, for the cases 
(I) and (II). The numerical errors demonstrate reasonable accuracy in both 
cases, especially in the case of using the scheme (II). 
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Table 1: Numerical errors produced by the new method using scheme (I). 


M N Maxerror RMSE M WN _ Maxerror RMSE 
10 10 3.37E-2 8.40E-3 30 15 9.81E-3 4.80E-3 
20 10 1.43E-2 6.89E-3 40 20 7.40E-3 3.70E-3 
30 =610 1.51E-2 7.68E-3 50 25 6.01E-3 2.95E-3 
40 10 1.54E-2 8.02E-3 60 30 4.91E-3 2.50E-3 
50 10 1.56E-2 8.20E-3 70 35 4.27E-3 2.22E-3 


© Analytical solution 
+ Present method 


1 1 1 L L L L 1 1 1 L L L L 1 1 1 L L 
1 08 06 O04 02 6 O2 O4 O68 O08 1 Ol —gs a6 a4 02 0 02 G4 08 08 1 
x x 


Figure 2: Comparing the exact solution and the numerical solution by scheme (I) with 
M = 70 and N = 35 and configuration of the error function. 


In order to show the efficiency of the new method, we compare the nu- 
merical results with those presented in [13]. The errors obtained by scheme 
(II) and those presented in [13] are listed in Table 3 using the same mesh 
size and time steps. It is observed that although in some cases our method 
produces more accurate solutions, overall, the difference between the errors 
is not significant. However, the new method has some advantages over the 
previous ones, that is, iterative methods. In the previous works, an initial so- 
lution should be carefully selected, otherwise the convergence of the iterative 
method is not guaranteed. In addition, as seen in Table 3, in the previous 
work, the solution is obtained at the expense of performing some iterations, 
which considerably increase the computational costs. Moreover, a parameter 
is needed in the iterative method (see [13]). Choosing this parameter is also 
important in order to achieve convergence. 


It should be noted, in the suggested method, a global system of equations 
is solved rather than the local ones. Although, apparently, this may increase 
the computational costs, the resulting system has a sparse coefficient matrix, 
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Table 2: Numerical errors produced by the new method using scheme (II). 


M  N_ Maxerror RMSE M WN _ Maxerror RMSE 
10 10 3.10E-2 8.53E-3 30 15 3.85E-3 9.12E-4 
20 10 7.48E-3 2.01E-3 40 20 2.36E-3 5.18E-4 
30 ©6110 3.50E-3 118E-3 50 25 1.62E-3 3.33E-4 
40 10 2.86E-3 1.10E-3 60 30 1.18E-3 2.32E-4 
50 10 2.61E-3 1.00E-3 70 35 8.75E-4 1.71E-4 
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which can be treated by the efficient solvers. Moreover, the iterations required 
in the case of using local matrices, are avoided. 
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1 1 1 L L 
1 08 06 O04 02 0 
x 


L L 1 1 
02 04 O68 o8 1 


1. 


L L L L 1 
8 06 04 02 0 


1 L L L 
02 04 06 o8 1 


Figure 3: Comparing the exact solution and numerical solution by scheme (II) with 


M = 70, N = 35 and configuration of the error function. 


6 Conclusions 


A nonoverlapping DDM was applied to the FBHE in a one-dimensional case. 
The work was based on partitioning the spatial domain and considering each 
part as an independent forward or a backward subproblem. A forward and a 
backward finite difference schemes were employed for local problems followed 
by assembling the equations in a global one. The under determined system 
of equations, arising from the lack of the interface boundary conditions, was 
treated by a physical matching condition and a square system was obtained. 
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Table 3: Comparing the error values (Max error) of scheme (II) with the previous results 
produced by the iterative method with K denoting the number of iterations. 


M N kK ~ Iterative method New method 


4 4 24 1.36E-2 7.06E-2 
8 16-338 3.68E-3 1.02E-2 
16 64 63 9.22E-3 1.10E-3 
32 256 100 2.42E-4 2.56E-4 


Dealing with a single global system of equations, while avoiding iterations 
used in the previous attempts, led to increasing the computational efficiency 
and keeping a reasonable accuracy for the solution. The proposed method 
performed well in the one-dimensional case. The application of the present 
method to the case of two-dimensional is also under consideration. 
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